suppressPackageStartupMessages({
    library(Seurat)
    library(cowplot)
    library(ggplot2)
    library(pheatmap)
    library(rafalib)
    library(clustree)
})
 
library(future)
plan("multiprocess", workers = 8)
options(future.globals.maxSize = 48000 * 1024^2)

1 Graph clustering

if(!file.exists("../analysis/filtered_IMM_DN_int_clus.rds")){
  alldata = readRDS(file = "../analysis/filtered_IMM_DN_int.rds")
  alldata@active.assay = "CCA"

alldata <- FindNeighbors(alldata, dims = 1:30, k.param = 350, prune.SNN = 1/15)
names(alldata@graphs)

pheatmap(alldata@graphs$CCA_nn[1:200, 1:200], col = c("white", "black"), border_color = "grey90", 
    legend = F, cluster_rows = F, cluster_cols = F, fontsize = 2)

# Clustering with louvain (algorithm 1)
for (res in c(0.1, 0.25, 0.5, 1, 1.5, 2)) {
    alldata <- FindClusters(alldata, graph.name = "CCA_snn",
                            resolution = res, algorithm = 1)
}


saveRDS(alldata, file = "../analysis/filtered_IMM_DN_int_clus.rds")

}else{
  alldata = readRDS(file = "../analysis/filtered_IMM_DN_int_clus.rds")
}
plot_grid(ncol = 2,
          DimPlot(alldata, reduction = "umap", group.by = "CCA_snn_res.0.1", 
                  label = TRUE) + NoLegend() + 
            ggtitle("louvain_0.1") + theme(legend.position = "bottom"), 
          DimPlot(alldata, reduction = "umap", group.by = "CCA_snn_res.0.25", 
                  label = TRUE) + NoLegend() + 
            ggtitle("louvain_0.25") + theme(legend.position = "bottom"), 
          DimPlot(alldata, reduction = "umap", group.by = "CCA_snn_res.0.5", 
                  label = TRUE) + NoLegend() + 
            ggtitle("louvain_0.5") + theme(legend.position = "bottom"),
          DimPlot(alldata, reduction = "umap", group.by = "CCA_snn_res.1", 
                  label = TRUE) + NoLegend() + 
            ggtitle("louvain_1") + theme(legend.position = "bottom"))

clustree(alldata@meta.data, prefix = "CCA_snn_res.")

1.1 QC, CCA louvain 0.25

table(alldata$CCA_snn_res.0.25)

VlnPlot(alldata, features = c("nFeature_CCA", "nCount_CCA", "percent_mito", "percent_ribo", 
                              "S.Score", "G2M.Score"), group.by = "CCA_snn_res.0.25" , 
    ncol = 3, pt.size = 0.1)

VlnPlot(alldata, features = c("nFeature_CCA", "nCount_CCA", "percent_mito", "percent_ribo", 
                              "S.Score", "G2M.Score"), group.by = "CCA_snn_res.0.25" , 
    ncol = 3, pt.size = 0)

#plot as proportion or percentage of cluster
ggplot(alldata@meta.data, aes(x=CCA_snn_res.0.25, fill=orig.ident)) + geom_bar(position = "fill")

#plot as proportion or percentage of cluster
ggplot(alldata@meta.data, aes(x=orig.ident, fill=CCA_snn_res.0.25)) + geom_bar(position = "fill")

ggplot(alldata@meta.data, aes(x=Treatment, fill=CCA_snn_res.0.25)) + geom_bar(position = "fill") + facet_wrap(~Type)

ggplot(alldata@meta.data, aes(fill=Type, x=CCA_snn_res.0.25)) + geom_bar(position = "fill")

DotPlot(alldata, features = c("Ptprc","Epcam"), group.by = "CCA_snn_res.0.25",
    assay = "CCA") + coord_flip()

## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 11520  7328  6802  5322  5179  4375  3126  2667  2616  2198  1918  1620  1546 
##    13    14    15    16 
##  1492  1312   567   289

1.2 QC, CCA louvain 0.5

table(alldata$CCA_snn_res.0.5)

VlnPlot(alldata, features = c("nFeature_CCA", "nCount_CCA", "percent_mito", "percent_ribo", 
                              "S.Score", "G2M.Score"), group.by = "CCA_snn_res.0.5" , 
    ncol = 3, pt.size = 0.1)

VlnPlot(alldata, features = c("nFeature_CCA", "nCount_CCA", "percent_mito", "percent_ribo", 
                              "S.Score", "G2M.Score"), group.by = "CCA_snn_res.0.5" , 
    ncol = 3, pt.size = 0)

#plot as proportion or percentage of cluster
ggplot(alldata@meta.data, aes(x=CCA_snn_res.0.5, fill=orig.ident)) + geom_bar(position = "fill")

#plot as proportion or percentage of cluster
ggplot(alldata@meta.data, aes(x=orig.ident, fill=CCA_snn_res.0.5)) + geom_bar(position = "fill")

ggplot(alldata@meta.data, aes(x=Treatment, fill=CCA_snn_res.0.5)) + geom_bar(position = "fill") + facet_wrap(~Type)

ggplot(alldata@meta.data, aes(fill=Type, x=CCA_snn_res.0.5)) + geom_bar(position = "fill")

DotPlot(alldata, features = c("Ptprc","Epcam"), group.by = "CCA_snn_res.0.5",
    assay = "CCA") + coord_flip()

## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
## 7323 7024 6791 5324 5179 4469 4375 3126 2667 2616 1920 1622 1563 1562 1498 1325 
##   16   17   18 
##  634  570  289

1.3 QC, CCA louvain 1

table(alldata$CCA_snn_res.1)

VlnPlot(alldata, features = c("nFeature_CCA", "nCount_CCA", "percent_mito", "percent_ribo", 
                              "S.Score", "G2M.Score"), group.by = "CCA_snn_res.1" , 
    ncol = 3, pt.size = 0.1)

VlnPlot(alldata, features = c("nFeature_CCA", "nCount_CCA", "percent_mito", "percent_ribo", 
                              "S.Score", "G2M.Score"), group.by = "CCA_snn_res.1" , 
    ncol = 3, pt.size = 0)

#plot as proportion or percentage of cluster
ggplot(alldata@meta.data, aes(x=CCA_snn_res.1, fill=orig.ident)) + geom_bar(position = "fill")

#plot as proportion or percentage of cluster
ggplot(alldata@meta.data, aes(x=orig.ident, fill=CCA_snn_res.1)) + geom_bar(position = "fill")

ggplot(alldata@meta.data, aes(x=Treatment, fill=CCA_snn_res.1)) + geom_bar(position = "fill") + facet_wrap(~Type)

ggplot(alldata@meta.data, aes(fill=Treatment, x=CCA_snn_res.1)) + geom_bar(position = "fill")

DotPlot(alldata, features = c("Ptprc","Epcam"), group.by = "CCA_snn_res.1",
    assay = "CCA") + coord_flip()

## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
## 6780 6745 6137 5325 5183 4480 4375 2670 2616 2411 1920 1627 1576 1563 1502 1341 
##   16   17   18   19   20   21 
##  857  704  637  570  569  289
sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS/LAPACK: /Users/jenfra/miniconda3/envs/singlecell/lib/libopenblasp-r0.3.15.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] future_1.27.0      clustree_0.5.0     ggraph_2.0.5       rafalib_1.0.0     
##  [5] pheatmap_1.0.12    ggplot2_3.3.6      cowplot_1.1.1      sp_1.5-0          
##  [9] SeuratObject_4.1.0 Seurat_4.1.1      
## 
## loaded via a namespace (and not attached):
##   [1] Rtsne_0.15            colorspace_2.0-3      deldir_1.0-6         
##   [4] ellipsis_0.3.2        ggridges_0.5.3        rstudioapi_0.13      
##   [7] spatstat.data_2.2-0   farver_2.1.0          leiden_0.3.9         
##  [10] listenv_0.8.0         graphlayouts_0.8.0    ggrepel_0.9.1        
##  [13] fansi_1.0.3           codetools_0.2-18      splines_4.0.5        
##  [16] cachem_1.0.6          knitr_1.37            polyclip_1.10-0      
##  [19] jsonlite_1.8.0        ica_1.0-2             cluster_2.1.2        
##  [22] png_0.1-7             rgeos_0.5-9           uwot_0.1.11          
##  [25] ggforce_0.3.3         shiny_1.7.2           sctransform_0.3.3    
##  [28] spatstat.sparse_2.1-1 compiler_4.0.5        httr_1.4.2           
##  [31] backports_1.4.1       assertthat_0.2.1      Matrix_1.4-1         
##  [34] fastmap_1.1.0         lazyeval_0.2.2        cli_3.3.0            
##  [37] tweenr_1.0.2          later_1.3.0           htmltools_0.5.3      
##  [40] tools_4.0.5           igraph_1.3.4          gtable_0.3.0         
##  [43] glue_1.6.2            RANN_2.6.1            reshape2_1.4.4       
##  [46] dplyr_1.0.9           Rcpp_1.0.9            scattermore_0.8      
##  [49] jquerylib_0.1.4       vctrs_0.4.1           nlme_3.1-155         
##  [52] progressr_0.10.0      lmtest_0.9-40         spatstat.random_2.2-0
##  [55] xfun_0.31             stringr_1.4.0         globals_0.16.0       
##  [58] mime_0.12             miniUI_0.1.1.1        lifecycle_1.0.1      
##  [61] irlba_2.3.5           goftest_1.2-3         MASS_7.3-55          
##  [64] zoo_1.8-9             scales_1.1.1          tidygraph_1.2.0      
##  [67] spatstat.core_2.4-4   promises_1.2.0.1      spatstat.utils_2.3-1 
##  [70] parallel_4.0.5        RColorBrewer_1.1-2    yaml_2.3.5           
##  [73] reticulate_1.25       pbapply_1.5-0         gridExtra_2.3        
##  [76] sass_0.4.0            rpart_4.1.16          stringi_1.7.8        
##  [79] highr_0.9             checkmate_2.1.0       rlang_1.0.4          
##  [82] pkgconfig_2.0.3       matrixStats_0.62.0    evaluate_0.15        
##  [85] lattice_0.20-45       tensor_1.5            ROCR_1.0-11          
##  [88] purrr_0.3.4           labeling_0.4.2        patchwork_1.1.1      
##  [91] htmlwidgets_1.5.4     tidyselect_1.1.2      parallelly_1.32.1    
##  [94] RcppAnnoy_0.0.19      plyr_1.8.7            magrittr_2.0.3       
##  [97] bookdown_0.24         R6_2.5.1              generics_0.1.3       
## [100] DBI_1.1.3             withr_2.5.0           mgcv_1.8-39          
## [103] pillar_1.8.0          fitdistrplus_1.1-8    survival_3.3-1       
## [106] abind_1.4-5           tibble_3.1.8          future.apply_1.8.1   
## [109] crayon_1.5.1          KernSmooth_2.23-20    utf8_1.2.2           
## [112] spatstat.geom_2.4-0   plotly_4.10.0         rmarkdown_2.13       
## [115] viridis_0.6.2         grid_4.0.5            data.table_1.14.2    
## [118] digest_0.6.29         xtable_1.8-4          tidyr_1.2.0          
## [121] httpuv_1.6.5          munsell_0.5.0         viridisLite_0.4.0    
## [124] bslib_0.4.0